#include "includes/basis-function.hh"

BasisFunction::BasisFunction(){
}

BasisFunction::~BasisFunction(){
}

void BasisFunction::getInitComponent(){
  for(int tmpi=0; tmpi<memory::phynumber;tmpi++){
    for(int tmpj=0; tmpj<memory::dgorder;tmpj++){
    }
  }
}

void BasisFunction::getDBasisLeft(double dbl[]){
  dbl[0]=0.0;
  dbl[1]=1.0;
  dbl[2]=-3.0;
  dbl[3]=6.0;
  dbl[4]=-10.0;      
  dbl[5]=15.0;      
  dbl[6]=-21.0;      
  dbl[7]=28.0;
}

void BasisFunction::getDBasisRight(double dbr[]){
  dbr[0]=0.0;
  dbr[1]=1.0;
  dbr[2]=3.0;
  dbr[3]=6.0;
  dbr[4]=10.0;      
  dbr[5]=15.0;      
  dbr[6]=21.0;      
  dbr[7]=28.0;      
}

void BasisFunction::getDBasisMiddle(double dbm[]){
  dbm[0]=0.0;
  dbm[1]=1.0;
  dbm[2]=0.0;
  dbm[3]=-3.0/2.0;
  dbm[4]=0.0;      
  dbm[5]=15.0/8.0;      
  dbm[6]=0.0;      
  dbm[7]=-35.0/16.0;      
}

void BasisFunction::getBasisLeft(double bl[]){
  bl[0]=1.0;
  bl[1]=-1.0;
  bl[2]=1.0;
  bl[3]=-1.0;
  bl[4]=1.0;
  bl[5]=-1.0;
  bl[6]=1.0;
  bl[7]=-1.0;  
}

void BasisFunction::getBasisRight(double br[]){
  br[0]=1.0;
  br[1]=1.0;
  br[2]=1.0;
  br[3]=1.0;
  br[4]=1.0;
  br[5]=1.0;
  br[6]=1.0;
  br[7]=1.0;  
}

void BasisFunction::getBasisMiddle(double bm[]){
  bm[0]=1.0;
  bm[1]=0.0;
  bm[2]=-1.0/2.0;
  bm[3]=0.0;
  bm[4]=3.0/8.0;
  bm[5]=0.0;
  bm[6]=-5.0/16.0;
  bm[7]=0.0;
}

void BasisFunction::getVarLeft(int i,double varl[]){
  for(int tmpi=0; tmpi<memory::phynumber;tmpi++)
    varl[tmpi]=0.0;
  for(int tmpi=0; tmpi<memory::phynumber;tmpi++){
    for(int tmpj=0; tmpj<memory::dgorder;tmpj++){
      varl[tmpi]+=pow(-1,memory::dgorder)*memory::variable[tmpj*memory::phynumber+tmpi][i];
    }
  }  
}

void BasisFunction::getVarRight(int i,double varr[]){
  for(int tmpi=0; tmpi<memory::phynumber;tmpi++)
    varr[tmpi]=0.0;
  for(int tmpi=0; tmpi<memory::phynumber;tmpi++){
    for(int tmpj=0; tmpj<memory::dgorder;tmpj++){
      varr[tmpi]+=memory::variable[tmpj*memory::phynumber+tmpi][i];
    }
  }  
}

void BasisFunction::getVarMiddle(int i,double varm[]){
  for(int tmpi=0; tmpi<memory::phynumber;tmpi++)
    varm[tmpi]=0.0;
  for(int tmpi=0; tmpi<memory::phynumber;tmpi++){
    varm[tmpi]+=memory::variable[0*memory::phynumber+tmpi][i];
    if (memory::dgorder>2 && memory::dgorder<=4)
      varm[tmpi]+=-1.0/2.0*memory::variable[2*memory::phynumber+tmpi][i];
    if (memory::dgorder>4 && memory::dgorder<=6)
      varm[tmpi]+=3.0/8.0*memory::variable[4*memory::phynumber+tmpi][i];
    if (memory::dgorder>6 && memory::dgorder<=8)
      varm[tmpi]+=-5.0/16.0*memory::variable[6*memory::phynumber+tmpi][i];
  }
}

void BasisFunction::getOnePhyLeft(int i,int phyi,double varl[]){
  varl[phyi]=0.0;
  // std::cout << "get one physical variable on the left" << std::endl;
  // std::cout << "phyi: " << phyi << std::endl;
  for(int tmpj=0; tmpj<memory::dgorder;tmpj++)
    varl[phyi]+=pow(-1,memory::dgorder)*memory::variable[tmpj*memory::phynumber+phyi][i];
  // varl[phyi]+=pow(-1,memory::dgorder)*memory::variable[phyi*memory::phynumber+tmpj][i];  
}

void BasisFunction::getOnePhyRight(int i,int phyi,double varr[]){
    varr[phyi]=0.0;
    // std::cout << "get one physical variable on the right" <<std::endl;
    for(int tmpj=0; tmpj<memory::dgorder;tmpj++)
      varr[phyi]+=memory::variable[tmpj*memory::phynumber+phyi][i];
    // varr[phyi]+=memory::variable[phyi*memory::phynumber+tmpj][i];    
}

void BasisFunction::getOnePhyMiddle(int i,int phyi,double varm[]){
  varm[phyi]=0.0;
  varm[phyi]+=memory::variable[0*memory::phynumber+phyi][i];
  if (memory::dgorder>2 && memory::dgorder<=4)
    varm[phyi]+=-1.0/2.0*memory::variable[2*memory::phynumber+phyi][i];
  if (memory::dgorder>4 && memory::dgorder<=6)
    varm[phyi]+=3.0/8.0*memory::variable[4*memory::phynumber+phyi][i];
  if (memory::dgorder>6 && memory::dgorder<=8)
    varm[phyi]+=-5.0/16.0*memory::variable[6*memory::phynumber+phyi][i];
}

// double BasisFunction::toxi(double x, double xcenter, double pillarlen){
//   return (x-xcenter)/pillarlen*2.0;
// }

void BasisFunction::getBasis(double binside[]){
  for(int tmpj=0; tmpj<memory::dgorder; tmpj++)
    for(int tmpk=0; tmpk<memory::intergnumber; tmpk++){
      // std::cout << -1.0+2.0*(tmpk+0.5)/memory::intergnumber << std::endl;
      // binside[tmpj*memory::intergnumber+tmpk]=legendreP(-1.0+2.0*(tmpk+0.5)/memory::intergnumber,tmpj);
      binside[tmpj*memory::intergnumber+tmpk]=legendreP(-1.0+2.0*tmpk/(memory::intergnumber-1.0),tmpj);
    }
}

//CHECK: Done-1
double BasisFunction::legendreP(double xi, int n){
  double res=0.0;
  if(n==0){
    res=1.0;
  } else if (n==1) {
    res=xi;
  } else if (n==2) {
    res=1.0/2.0*(-1.0+3.0*xi*xi);
  } else if (n==3) {
    res=1.0/2.0*(-3.0*xi+5.0*xi*xi*xi);
  } else if (n==4) {
    res=1.0/8.0*(3.0-30.0*xi*xi+35.0*pow(xi,4));
  } else if (n==5) {
    res=1.0/8.0*(15.0*xi-70.0*xi*xi*xi+63.0*pow(xi,5));
  } else if (n==6) {
    res=1.0/16.0*(-5.0+105.0*xi*xi-315.0*pow(xi,4)+231.0*pow(xi,6));
  }
  return res;
}

// CHECK: Done-1
double BasisFunction::dlegendreP(double xi, int n){
  double res=0.0;
  if(n==0){
    res=0.0;
  } else if (n==1) {
    res=1.0;
  } else if (n==2) {
    res=3.0*xi;
  } else if (n==3) {
    res=1.0/2.0*(-3.0+15.0*xi*xi);
  } else if (n==4) {
    res=1.0/8.0*(-60.0+140.0*xi*xi*xi);
  } else if (n==5) {
    res=1.0/8.0*(15.0-210.0*xi*xi+315.0*pow(xi,4));
  } else if (n==6) {
    res=1.0/16.0*(210.0-1260.0*xi*xi*xi+1386.0*pow(xi,5));
  }
  return res;
}

BasisFunction basis;
